
# Set up the packages, options, and working directory path
pacman::p_load(tidyverse, haven, RColorBrewer, binsreg)
set.seed(42)
options(scipen=999)

if (Sys.info()[['user']] == "alal") {
  root = "/home/alal/Dropbox/1_Research/ElecAdminFunding/ElecAdminFundingReplication"
  theme_set(lal_plot_theme())
} else {
  root = "~/Dropbox/ElecAdminFunding/ElecAdminFundingReplication"
  theme_set(theme_minimal())
}

blues = brewer.pal(8, 'Blues')[c(5, 8)]
partisan_colors = c(brewer.pal(8, 'Reds')[8],
  brewer.pal(8, 'Blues')[8])


###
# Set up the plot-making functions
###

# Make a function that prepares data for a binned scatter plot
bin_data = function(data, y, x){

  # Set up a simple dataset with the x and y variables
  dat = data %>%
    select(all_of(c(y, x))) %>%
    rename("y"=!!y, "x"=!!x)

  # Let binsreg choose the number of bins
  bins = binsregselect(dat$y, dat$x)$nbinsrot.regul
  n = nrow(dat)

  # Make the simple binscatter
  binned_dat = dat %>%
    # Divide the data into equal-sized bins based on the
    # auto-selected number of bins
    arrange(x) %>%
    mutate(bin = ceiling(row_number()/ceiling(n()/bins))) %>%
    # Get the bin means
    group_by(bin) %>%
    mutate(y_mean = mean(y),
      x_mean = mean(x),
      y_mean_first = ifelse(row_number()==1, y_mean, NA),
      x_mean_first = ifelse(row_number()==1, x_mean, NA)) %>%
    ungroup()

  return(list(binned_dat=binned_dat, bins=bins, n=n))
}


# Make a binscatter function
binscatter = function(data, y, x, xlab = "", ylab = "",
  xmin = NULL, xmax = NULL, ymin=NULL, ymax=NULL, xtext, ytext){

  #
  bin_output = bin_data(data=data, y=y, x=x)
  n_text = paste0(format(bin_output$n, nsmall=1, big.mark=","), " Counties")
  bin_size_text = paste0(round(bin_output$n/bin_output$bins), " Counties per Bin")

  # Make the simple binned scatter plot
  fig = bin_output$binned_dat %>%
    ggplot() +
      geom_point(aes(x=x_mean_first, y=y_mean_first), size=5, color="#08519C", alpha=0.5) + # #3182BD #6BAED6
      geom_smooth(aes(x=x, y=y), linewidth=1.5, color="#08519C", method="lm", se=FALSE) +
      labs(x=xlab, y=ylab) +
      theme(text = element_text(size=20, face="bold", color="black")) +
      theme(panel.grid.minor = element_blank()) +
      annotate("text", x=xtext, y=ytext + 5, label=n_text, size=6) +
      annotate("text", x=xtext, y=ytext, label=bin_size_text, size=6)
  if(!is.null(xmin)) fig  = fig + xlim(xmin, xmax)
  if(!is.null(ymin)) fig  = fig + ylim(ymin, ymax)
  return(fig)
}


###
# Set up the data
###

# Set up the analysis file
uvbm_states = c("WA", "CO", "OR", "HI", "UT")
ctcl = haven::read_dta(file.path(root, "/modified_data/ctcl.dta")) %>%
  mutate(treat_pct = treat*100,
    pop_terc = ceiling(percent_rank(vap)*3),
    pop_terc = ifelse(pop_terc==0, 1, pop_terc),
    vbm_change = ifelse(state %in% uvbm_states, 0, no_excuse))

# Create a file of trends on turnout and Dem vote share
trends = ctcl %>%
  group_by(treat, policy_year) %>%
  summarize(dem_vs = mean(vs_dem_pres), 
    turnout = mean(turnout_pres), .groups="drop") %>%
  mutate(treat = ifelse(treat==1, "Received Grant", 
    "Did Not Receive Grant"))


###
# Make the simple binscatter selection plots
###

# Make the 2016 dem vote share selection plot
fig_vs_selection = ctcl %>%
  filter(policy_year==2016, !is.na(treat),
    !is.na(vs_dem_pres)) %>%
  binscatter(y="treat_pct", x="vs_dem_pres",
    xlab = "Democratic Vote Share (2016)",
    ylab = "Applied for and Received Grant (%)",
    xmin=0, xmax=100, ymin=0, ymax=100,
    xtext=87, ytext=4)
ggsave(file.path(root, "output/vs_selection.pdf"), plot=fig_vs_selection,
  width = 10, height = 8)


###
# Make the simple trends plots
###

# Make the Dem vote share trends plot
trends_dem_vs = trends %>%
  ggplot(aes(x = policy_year, y = dem_vs, group = treat, 
      colour = as.factor(treat))) +
    geom_line(linewidth=1.5) + geom_point(size=2) +
    scale_color_manual(values = blues, name = "") +
    scale_x_continuous(limits = c(1991, 2021),
          breaks = seq(1992, 2020, by=4)) +
    geom_vline(xintercept = 2018, linetype = 'dotted') +
    labs(y = "Dem Vote Share (%)", x = "Year", colour = "Applied for and Received Grant") +
    theme(text = element_text(size=20, face="bold", color="black"),
          legend.position="bottom", legend.title= element_blank(),
          panel.grid.minor = element_blank())
ggsave(file.path(root, "output/trends_dem_vs.pdf"),
  plot = trends_dem_vs, height = 8, width = 10)

# Make the turnout trends plot
trends_turnout = trends %>%
  ggplot(aes(x = policy_year, y = turnout, group = treat, 
      colour = as.factor(treat))) +
    geom_line(linewidth=1.5) + geom_point(size=2) +
    scale_color_manual(values = blues, name = "") +
    scale_x_continuous(limits = c(1991, 2021),
          breaks = seq(1992, 2020, by=4)) +
    geom_vline(xintercept = 2018, linetype = 'dotted') +
    labs(y = "Turnout (%)", x = "Year", colour = "Received Grant") +
    theme(text = element_text(size=20, face="bold", color="black"),
          legend.position="bottom", legend.title= element_blank(),
          panel.grid.minor = element_blank())
ggsave(file.path(root, "output/trends_turnout.pdf"),
  plot = trends_turnout, height = 8, width = 10)


###
# Plot grant receipt by state
###

# Make the treatment by state bar plot
treat_by_state = ctcl %>%
  filter(policy_year==2020) %>%
  group_by(state) %>%
  summarize(treat_share = mean(100*treat)) %>%
  ggplot(aes(y = reorder(state, treat_share), x = treat_share)) +
    geom_col(color=blues[2], fill=blues[2]) +
    labs(x = "Share Receiving Grant (%)", y="") +
    theme(text = element_text(size=12, face="bold", color="black"))
ggsave(file.path(root, "output/treat_by_state.pdf"), treat_by_state, # device = cairo_pdf, 
  width = 8, height = 8)

